/* ====================================================
#   Copyright (C) 2021  All rights reserved.
#
#   Author        : Jinjing Li
#   Email         : josiahleeoaa@outlook.com
#   File Name     : libSpFission.C
#   Last Modified : 2025-05-18 17:07
#   Describe      :
#
# ====================================================*/

#include <cstdio>
#include <cmath>
#include "libSpFission.h"

using namespace std;
using namespace CLHEP;

void generator_spfission_basedon_freya_interface(int isotope, unsigned int nEvents, FILE *stream)
{
    double time = 0;
    int corr_mode = 3;
#ifdef OUTPUT_2_ROOT_FILE
    TFile *out_file = new TFile("output.root", "RECREATE");
    TH1D *h_nnu = new TH1D("h_nnu", "Neutron number", max_fold, -0.5, max_fold - 0.5);
    TH1D *h_pnu = new TH1D("h_pnu", "Gamma number", max_fold, -0.5, max_fold - 0.5);
    TH1D *h_neng[max_fold];
    TH1D *h_peng[max_fold];
    for (int i = 0; i < max_fold; i++)
    {
        h_neng[i] = new TH1D(Form("h_neng_%d", i + 1), "Neutron energy", 1000, 0, 20);
        h_peng[i] = new TH1D(Form("h_peng_%d", i + 1), "Gamma energy", 1000, 0, 20);
    }
#endif

    // set correlation between neutron and gamma: number and energy
    setcorrel_(&corr_mode);

    // start by printing some information to comment lines
    fprintf(stream, "# File generated by u_generator made by Jinjing Li.\n");

    for (unsigned int i = 0; i < nEvents; i++)
    {
        genspfissevt_(&isotope, &time);
        int n_neutron = getnnu_();
        int n_gamma = getpnu_();

        fprintf(stream, "%d\n", n_neutron + n_gamma); // number of particles from this vertex

        for (int j = 0; j < n_neutron; j++)
        {
            double eng_neutron = getneng_(&j);
            double age_neutron = getnage_(&j);
            double vecs[3] = {getndircosu_(&j), getndircosv_(&j), getndircosw_(&j)};
            double momentum = sqrt(pow(eng_neutron / GeV + neutron_mass_c2 / GeV, 2) - pow(neutron_mass_c2 / GeV, 2));
            dirs this_dir(vecs);

#ifdef OUTPUT_2_ROOT_FILE
            h_neng[j]->Fill(eng_neutron);
#endif

            fprintf(stream, "1\t2112 0 0 %e %e %e %e %e\n", momentum * (this_dir.x_norm),
                    momentum * (this_dir.y_norm),
                    momentum * (this_dir.z_norm),
                    neutron_mass_c2 / GeV, age_neutron);
        }
        for (int j = 0; j < n_gamma; j++)
        {
            double eng_gamma = getpeng_(&j);
            double age_gamma = getpage_(&j);
            double vecs[3] = {getpdircosu_(&j), getpdircosv_(&j), getpdircosw_(&j)};
            double momentum = eng_gamma / GeV;
            dirs this_dir(vecs);

#ifdef OUTPUT_2_ROOT_FILE
            h_peng[j]->Fill(eng_gamma);
#endif

            fprintf(stream, "1\t22 0 0 %e %e %e %e %e\n", momentum * (this_dir.x_norm),
                    momentum * (this_dir.y_norm),
                    momentum * (this_dir.z_norm),
                    0., age_gamma);
        }
#ifdef OUTPUT_2_ROOT_FILE
        h_nnu->Fill(n_neutron);
        h_pnu->Fill(n_gamma);
#endif
    }

#ifdef OUTPUT_2_ROOT_FILE
    h_nnu->Write();
    h_pnu->Write();
    for (int i = 0; i < max_fold; i++)
    {
        h_neng[i]->Write();
        h_peng[i]->Write();
    }
    out_file->Close();
#endif
}

void generator_spfission_basedon_freya_fissevent(int isotope, unsigned int nEvents, FILE *stream)
{
    double energy_MeV = 0;
    double nubar = 2.01;
    double time = 0.;

    char errors[MAX_ERROR_LENGTH];

#ifdef OUTPUT_2_ROOT_FILE
    TFile *out_file = new TFile("output_fissionevent.root", "RECREATE");
    TH1D *h_nnu = new TH1D("h_nnu", "Neutron number", max_fold, -0.5, max_fold - 0.5);
    TH1D *h_pnu = new TH1D("h_pnu", "Gamma number", max_fold, -0.5, max_fold - 0.5);
    TH1D *h_neng[max_fold];
    TH1D *h_peng[max_fold];
    for (int i = 0; i < max_fold; i++)
    {
        h_neng[i] = new TH1D(Form("h_neng_%d", i + 1), "Neutron energy", 1000, 0, 20);
        h_peng[i] = new TH1D(Form("h_peng_%d", i + 1), "Gamma energy", 1000, 0, 20);
    }
#endif

    init();

    // start by printing some information to comment lines
    fprintf(stream, "# File generated by u_generator by Jinjing Li.\n");

    for (unsigned int i = 0; i < nEvents; i++)
    {
        fissionEvent *fe = new fissionEvent(isotope, time, nubar, energy_MeV, 0);
        int errorlength = MAX_ERROR_LENGTH;
        fe->getFREYAerrors(&errorlength, &errors[0]);
        /*
        if (errorlength > 1)
        {
          printf("%s\n", errors);
          exit(1);
        }
        */

        int n_neutron = fe->getNeutronNu();
        int n_gamma = fe->getPhotonNu();
        fprintf(stream, "%d\n", n_neutron + n_gamma); // number of particles from this vertex

        // neutron
        for (int j = 0; j < n_neutron; j++)
        {
            double eng_neutron = fe->getNeutronEnergy(j);
            // you need to know: u2+v2+w2=1
            double u_neutron = fe->getNeutronDircosu(j);
            double v_neutron = fe->getNeutronDircosv(j);
            double w_neutron = fe->getNeutronDircosw(j);
            double age_neutron = fe->getNeutronAge(j);
            double momentum = sqrt(pow(eng_neutron / GeV + neutron_mass_c2 / GeV, 2) - pow(neutron_mass_c2 / GeV, 2));

#ifdef OUTPUT_2_ROOT_FILE
            h_neng[j]->Fill(eng_neutron);
#endif

            fprintf(stream, "1\t2112 0 0 %e %e %e %e %e\n", momentum * u_neutron, momentum * v_neutron, momentum * w_neutron, neutron_mass_c2 / GeV, age_neutron);
        }

        // gamma
        for (int j = 0; j < n_gamma; j++)
        {
            double eng_gamma = fe->getPhotonEnergy(j);
            // you need to know: u2+v2+w2=1
            double u_gamma = fe->getPhotonDircosu(j);
            double v_gamma = fe->getPhotonDircosv(j);
            double w_gamma = fe->getPhotonDircosw(j);
            double age_gamma = fe->getPhotonAge(j);
            double momentum = eng_gamma / GeV;
#ifdef OUTPUT_2_ROOT_FILE
            h_peng[j]->Fill(eng_gamma);
#endif
            fprintf(stream, "1\t22 0 0 %e %e %e %e %e\n", momentum * u_gamma, momentum * v_gamma, momentum * w_gamma, 0., age_gamma);
        }
#ifdef OUTPUT_2_ROOT_FILE
        h_nnu->Fill(n_neutron);
        h_pnu->Fill(n_gamma);
#endif

        delete fe;
    }
#ifdef OUTPUT_2_ROOT_FILE
    h_nnu->Write();
    h_pnu->Write();
    for (int i = 0; i < max_fold; i++)
    {
        h_neng[i]->Write();
        h_peng[i]->Write();
    }
    out_file->Close();
#endif
}

void init(void)
{
    unsigned short int s[3] = {1234, 5678, 9012};
    TRandom3 gr(0);
    s[0] = int(gr.Rndm() * 1e5);
    s[1] = int(gr.Rndm() * 1e5);
    s[2] = int(gr.Rndm() * 1e5);
    seed48(s);
    fissionEvent::setCorrelationOption(3);
    return;
}
